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Abstract 

We discuss the salient features of Zolotarev optimal rational ap- 
proximation for the inverse square root function, in particular, for its 
applications in lattice QCD with overlap Dirac quark. The theoretical 
error bound for the matrix-vector multiplication H^H^j^^Y is de- 
rived. We check that the error bound is always satisfied amply, for any 
QCD gauge configurations we have tested. An empirical formula for the 
error bound is determined, together with its numerical values ( by eval- 
uating elliptic functions ) listed in Table 2 as well as plotted in Figure 
3. Our results suggest that with Zolotarev approximation to (H^)^ 1 ^ 2 , 
one can practically preserve the exact chiral symmetry of the overlap 
Dirac operator to very high precision, for any gauge configurations on 
a finite lattice. 
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1 Introduction 



In lattice QCD with overlap Dirac quarks [1, 2], one encounters the challenging 
problem of taking the inverse square root of a positive definite Hermitian 
matrix, which stems from the overlap Dirac operator for massless fermion 



where H w denotes the Hermitian Wilson-Dirac operator with a negative pa- 
rameter — mo, 

H w = ^D W =75(-mo + 7A + W , (2) 

7^t M the naive fermion operator, and W the Wilson term. It is well-known 
that (1) is the simplest and the most elegant realization of chiral symmetry on 
a finite lattice [3] 

D l5 + l5 D = am^D^D , (3) 

and it also fulfils all basic requirements ( i.e., doublers-free, 75-hermiticity, 
exponentially-local, and correct continuum behavior ) for a decent lattice Dirac 
operator. However, in practice, how to compute physical quantities in lattice 
QCD with overlap Dirac operator (1) is still a challenging problem, due to the 
inverse square root of in (1), which cannot be evaluated analytically in a 
closed form, nor numerically via diagonalization since the memory requirement 
exceeds the physical memory of the present generation of computers. Thus, at 
this stage, it is necessary to replace (H 2 )' 1 / 2 with a good approximation, in 
any lattice QCD calculations with overlap Dirac quark. The relevant question 
is what is the optimal rational approximation r n (x) ( of degree n ) for the 
inverse square root function x~ 1//2 ,x G [1,6], in view of the convergence of a 
rational approximation can attain 



max 

l<x<b 



|1- y/Zr n (x)\ = Cl e- c * n , ci, c 2 >0, (4) 



which is much faster than that of any polynomial approximation. It turns out 
that the optimal solution has been obtained by Zolotarev [4] more than 100 
years ago. A detailed discussion of Zolotarev's result can be found in Akhiezer's 
two books [5, 6]. Unfortunately, Zolotarev's optimal rational approximation 
has been overlooked by the numerical algebra community until recent years 1 . 

A comparative study of Zolotarev's approximation versus other schemes for 
computing the product of the matrix sign function H^H"^)^ 1 / 2 with a vector 
Y, has been reported in Ref. [8]. 

1 See, for example, Ref. [7]. 
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Recently, we have used Zolotarev optimal rational approximation to com- 
pute overlap Dirac quark propagator in the quenched approximation [9]. The 
parameters of the pseudoscalar meson mass formula to one-loop order in chiral 
perturbation theory are determined, and from which the light quark masses 
m Ut d and m s can be extracted [10] with the experimental inputs of pion and 
kaon masses, and the pion decay constant. 

Nevertheless, in Refs. [8, 9], the underlying principles and salient features 2 
of Zolotarev optimal rational approximation have not been addressed. Also, 
in Refs. [7, 8, 9], only one option (70) of the two possible choices [(29) and 
(70)] of Zolotarev's rational polynomials has been exploited, however, these 
two options are complementary to each other, and in some cases it is more 
effective to use (29) than (70), in computing (H^)' 1 / 2 times a column vector 
Y, as discussed in Section 5. 

In this paper, we discuss the salient features of Zolotarev optimal rational 
approximation, in particular, for its application in lattice QCD with overlap 
Dirac operator. We start with the basic principles of rational approximation, 
and show that Zolotarev's rational polynomial is indeed the optimal rational 
approximation for the inverse square root function. Then we examine how 
well the Zolotarev optimal rational approximation (48) for the overlap Dirac 
operator can preserve the exact chiral symmetry (3) on a finite lattice. For 
the overlap Dirac operator, (3) is equivalent to 



Thus the question is how much the exact relation (5) is violated if one replaces 
(H^)" 1 / 2 with Zolotarev's optimal rational polynomial r^(if^). The chiral 
symmetry breaking can be measured in terms of the deviation 



It turns out that the deviation (6) has a theoretical upper bound (62), 
which is a function of n ( the degree of the rational polynomial ), and b = 
Knax I 'tfnim where ^max an d ^min denote the maximum and the minimum of 
the eigenvalues of H^. Thus, for any given gauge configuration, one can use 
the theoretical upper bound to determine what values of n and b ( i.e., how 
many low-lying eigenmodes of H 2 should be projected out ) are required to 
attain one's desired accuracy in preserving the chiral symmetry of the overlap 

2 In particular, the theoretical error bound for the product i?i„(.ff 2 ) _1 / 2 Y can be de- 
rived in terms of the deviation (35), as shown in Section 5. Moreover, we clarify that the 
coefficients do (49) and Dq (74) can be obtained explicitly, without using the condition : 
max[l - v / ^ r n( x )]li< ; E<b = -min[l - v/^n(a0]|i<x<6- 
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(6) 



yty 
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Dirac operator. In practice, one has no difficulties to achieve A z < 1CT 12 for 
any gauge configurations on a finite lattice. 

The outline of this paper is as follows. In Sections 2-4, we review the basic 
principles of rational approximation, and show that Zolotarev's rational poly- 
nomial is indeed the optimal rational approximation for the inverse square root 
function. In Section 5, we discuss the Zolotarev approximation for (i/ 2 ) -1 / 2 
in the overlap Dirac operator, and derive the theoretical error bound for the 
matrix- vector multiplication H^H^^^Y ( Y : any nonzero column vector 
), which is the most essential operation in computing the propagator of over- 
lap Dirac quark by nested conjugate gradient. We check that the theoretical 
error bound is always satisfied amply, for any QCD gauge configurations we 
have tested. Then the numerical values of the theoretical error bound are 
computed, and listed in Table 2, as well as plotted in Figure 3. An empirical 
formula for the error bound is determined from the numerical data. In Section 
6, we conclude with some remarks. 



2 de la Vallee-Poussin's theorem 

First we consider the following problem. Given any two positive and continuous 
functions f(x) and g(x) for x G [1,6], the problem is to find the irreducible 
rational polynomial of degree n with positive real coefficients, 

/ \ P{x) p n X n + Pn-lX 11 ' 1 H h p 

r(x) = — — = , (7) 

q(x) q n x n + q n -ix n - 1 + ■ ■ ■ + q 

such that the deviation of g(x)r(x) from f(x) is the minimum. Here the 
deviation is defined as the maximum of \f(x) — g(x)r(x)\ for the entire interval 
[1,6], i.e., 

d(r) = max \f(x) — g(x)r(x)\ . (8) 

x<=[l,b] 

Now if there exists r(x) such that the difference 

S( x ) = f(x) - g(x)r(x) (9) 

changes sign alternatively 2n + 2 times within the interval [1,6], and attains 
its maxima and minima, say, 

S(x) = Si, -5 2 , ■ ■ ■ , +52n+i, ~<W2 (Si>0) (10) 
at consecutive points, 

Xi < X 2 < ■ ■ ■ < X 2n +2 , (11) 
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then it can be shown that 

d(r) > min{£i, 5 2 , ■ • • , 5 2n+2 } = d min (12) 

for all irreducible rational polynomial r(x) as defined by (7). This can be 
asserted as follows. 

The strategy is to assume the contrary is true, and then show that it leads 
to contradiction. Suppose there exists an irreducible rational polynomial R(x) 
of degree n, such that d(R) < d min . Then 

D(x) = g(x)R(x) — g(x)r(x) 

= f(x) ~ g(x)r(x) - [f(x) - g(x)R(x)] 

= 5(x)-A(x) (13) 

is a continuous function, and it also changes sign alternatively at least 2n + 2 
times in the interval [1, b], since |A(x)| < d(R) < d m i n . Thus D(x) has at least 
2n + 1 zeros in the interval (1,6). On the other hand, 

D{x)=g(x)[R(x)-r(x)]=g{x)^ (14) 

where u(x) and v(x) are polynomials of degree 2n. Therefore D(x) cannot 
have more than 2n zeros, i.e., a contradiction. This completes the proof of de 
la Vallee-Poussin's theorem. 

In general, de la Vallee-Poussin's theorem[5] asserts that if there exists an 
irreducible rational polynomial of the form 

r M (x) = ^l^ 1 ^:''^ ,(™>n, > 0) (15) 
q m x m + q m ^ix m 1 + • • • + g 

such that 8{x) = f(x) — g(x)r^ n ' m \x) has n + m + 2 alternate change of sign 
in the interval [1,6], and attains its maxima and minima, say, 

5( x ) = 6 U -5 2 , • • • , (-l) n+m+1 5 n+m+2 (5 t > 0) 

at consecutive points, 

X\ < x 2 < ■ ■ ■ < x n+m + 2 , 

then 

d(r) > min{(5i, 5 2 , ■ ■ ■ , 5 n+m+2 } = d min 

for all irreducible rational polynomial r^ n ' m \x). 

In the next section, we show that there exists an optimal r z (x) such that 
d(rz) is minimum, i.e., Zolotarev's solution [4]. 
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3 Zolotarev's optimal rational approximation 

First, we recall some well-known properties of Jacobian elliptic functions [6]. 
The Jacobian elliptic function sn(w; k) — rj is defined by the elliptic integral 

u (rj) = / . . 16 

JO ^(l-t2)(l_ K 2 t 2) 

It is a meromorphic function of u and k, which is defined over C 2 . For r) = 1, 
(16) becomes 

r l dt 

«(1) = K — , (17) 

7 ° ^/(l-f 2 )(l-K 2 f 2 ) 

the complete elliptic integral of the first kind with modulus «. 

The important feature of sn(u; k) is that it is periodic in both real and 
imaginary parts of u, 

sn(u + AK; k) = sn(w; k) (18) 
sn(w + 2iK'; re) = sn(u; re) (19) 

where K' is complete elliptic integral of the first kind with modulus re' = 
Now if we define 

x(u; re) = sn 2 (w; re) , (20) 

then x has periods 2K and 2iK', since sn(w + 2K; re) = — sn(u; re). 

The crucial formula for Zolotarev's optimal rational approximation is the 
second principal n-th degree transformation of Jacobian elliptic function [6]. 
In terms of x (20), it reads 



/ ( u , \ l~, r 1 -A- 1 + x(u: k)/c2i ,^„ s 

' x tt; a ) = v^(«;«)t7TT; — / // 21 

VM / v v Mf} x 1 + x(u; «)/c2i_i 



where 



2n+1 e 2 (|^ ;«') 
A = n 7W(S)^ j A , (22) 



2 ( 



n sn 2 f (2LML. K >) 
M = n V , (23) 



and 6 denotes the elliptic theta function. Note that x(u/M; A) has periods 
2L = 2K/M and 2W = 2iK'/[(2n + 1)M], 

/ u 2iK' \ (u \ 

x \T7+ ,n„ I = * TF5 A • (26) 



M (2n+l)M' / VM 

Now restricting u = K + iv, we have 

, . n / \ 1 — sn 2 (?t>;re) 1 

x(u; re) = sn (u; re) = 3-—- = 7I — T , 27 

1 — re^sm(2t> ; re) 1 — re' sm(t>; re') 

where the identity 

™ 2 ^ = ~T^k < 28 > 

has been used. Equation (27) implies that x increases from 1 to 1/re 2 as 
v increases from to K', since sn(0; re') = 0, and sn(i^'; re') = 1, from the 
definition (16). 
Now consider 

2A 1 " l+x/c 2l 

rz(x) = TTxmE^J^ m 

as a rational approximation to 1/y/x. Then the deviation can be obtained 
from 

6 z (x) = 1 - Vtrz(x) = 1 - ^sn (*±™; a) (30) 

where (21) has been used. Using identity (27) with substitutions re <— A and 
u <— (K + iv)/M, then (30) becomes 

Sz(x) = 1 - 2A 1 , (31) 

1+ V i - A,28n2 (^ A ') 



where A' = y/l — A 2 . As v changes from to K', v/M changes from to 
K'/M = (2n + 1)V . This implies that sn 2 (v/M; A') attains its minima and 
maxima alternatively as 0,1, •••,0,1 at v/M = 0, L', ■ ■ • , 2nL', (2n + 1)1/. 
Correspondingly, 5z(x) reaches its maxima and minima alternatively as 

1_A 1-A _1^A 

1 + A' l + A''"'l + A' 1 + A 1 ' 
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at v/M = 0, L', • • • , 2nL', (2n + 1)L', which correspond to 2n + 2 successive 



That is, 



(33) 



1 - /t' 2 sn 2 ( i 1 ^ r ; K J 
From (32), we conclude that 

d(r z ) = = d min , (35) 

and Zolotarev's solution (29) is the optimal rational approximation for the 
inverse square root function, according to de la Vallee-Poussin's theorem. 

At this point, it is instructive to plot Sz(x) = 1 — y/xrz{x) explicitly, and 
to see how it attains its maxima and minima alternatively. This is shown in 
Fig. 1 for n = 6 and b = k~ 2 = 1000. Note that there are exactly 2n + 2 = 14 
alternate change of sign in [1, 1000], with maxima and minima at 

x = 1, 1.145, 1.664, 2.858, 5.415, 10.80, 22.05, 45.34, 
92.59, 184.7, 349.9, 600.9, 873.3, 1000. 

In Fig. 2, we plot the positions Xi (34) of the maxima and minima of 
Sz(x) = l—y/xrz(x) for n = 20 and b = hT 2 = 6000. Note that the distribution 
of Xi as shown in Fig. 2 is quite generic for any values of n and b. The maxima 
and minima near both ends ( x\ — 1 and rc2„ + 2 = b ) are densely packed even 
in the logarithm scale, while those at the central region seem to be evenly 
distributed in log scale ( i.e., exponentially in the linear scale ). 



4 Chebycheff's theorem 

Even though Zolotarev's rational approximation is optimal among all irre- 
ducible rational polynomials of degree n, with S(x) = f(x) — g(x)r(x) having 
2n + 2 alternate change of sign in [1,6], one may still ask whether there exists 
an optimal irreducible rational polynomial of the same degree but its S(x) has 
smaller number of alternate change of sign in [1,6]. Such a possibility is ruled 
out by Chebycheff's theorem, which can be shown as follows. 

Again, our strategy is to assume the contrary is true, and then show that 
it leads to contradiction. Suppose there exists an optimal rational polynomial 
R(x) = P(x)/Q(x) of degree n with deviation d(R) = d min , but its A(x) = 
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f(x) — g(x)R(x) only has m = 2n + 1 alternate change of sign in [1,6]. Now 
assume A(x) attains its maxima and minima, say, 

A(x) = Ai,-A 2 ,+A 3 ,.-.,+A m (A,>0) (36) 

at consecutive points 3 

1 = xi < x 2 < ■ ■ ■ < x m = b . (37) 

Then the interval [1, b] can be divided into m subintervals, 

[&,&]»••• ,[£m-i,b] , Xi<^<a: i+ i (38) 

such that the following two inequalities can be satisfied alternatively, 

-dmin + e < A(x) < rf min , x G [6i-2, i = 1, • • • , n + 1 (39) 
-d m in < A(x) < d min - e , x G fw] , i = 1, • • • , n (40) 

where e is a positive number which is much less than d min , and the notations 
£o — 1 an d £,m = b have been used. 

Next consider the following polynomial of degree In = m — 1 

S(x) = (x-^)(x-^)---(x-U~i) (41) 

which has sign change at £j, i = 1, • • • , m — 1. Using the fact that P(x) and 
Q(x) have no common factor, one can find polynomials U(x) and V(x) of 
degree n, with positive real coefficients, such that 

S(x) = Q(x)U(x) - P(x)V(x) . (42) 

Now consider the following irreducible rational polynomial of degree n 

w(x) = m±^m (43) 

W[X) Q(x) + aV(x) [A6) 
where a is a positive real parameter. Then, the difference 

f(x) - g(x)W(x) = f(x) - g(x)R(x) + g(x)R(x) - g(x)W(x) 

= A(x) - ag(x) Q( X )lQwiaV( X )] ' (44) 

is a function of a which can be chosen to be some tiny positive numbers such 
that the magnitude of the second term on the r.h.s. of (44) is much less than 
dmin - e for all x G [1,6], i.e., 

\S(x)\ 

° < a9{x) Q(x)[Q(x) + aV(x)] * » d ™ - " ■ (45) 



3 Without loss of generality, we set xi = 1, and x m = b. 
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Since A(x) and S(x) both have m = 2n + 1 alternate change of sign in 
[1,6], and the magnitude of the second term on r.h.s. of (44) is very small 
comparing to d min — e, it follows that f(x) — g(x)W(x) also has m = 2n + 1 
alternate change of sign in [1,6]. 

Then, using (39), (41), (44) and (45), we obtain 

-d m in + €-/!< f(x) - g(x)W (x) < d min - fj, , (46) 

for x G [C,2i-2, ^2j-i], i = 1, • ■ • ,n + 1. Similarly, we have 

-d m in + fi<f(x)- g(x)W(x) < d min -e + fi , (47) 

for x E [£2i-i, 6i], i = 1, ••• ,n . Thus we obtain that d(W) = d m in — A* < d min , 
a contradiction to our assumption that R(x) is the optimal one. This completes 
the proof of Chebycheff's theorem. 

In general, Chebycheff's theorem asserts that the optimal irreducible ratio- 
nal polynomial of the form (15) must satisfy the criterion that the difference 
S(x) = f(x) — g{x)r^ n,m \x) has n + m + 2 alternate change of sign in the 
interval [1,6]. 



5 Zolotarev approximation for the overlap 

In this section, we first derive the theoretical error bound for the the matrix- 
vector multiplication if w (.ff,^) -1 / 2 Y ( Y : any nonzero column vector ) with 
Zolotarev approximation for (if 2 ) -1 / 2 . Then the numerical values of the error 
bound are computed, and listed in Table 1, as well as plotted in Figure 3. An 
empirical formula for the error bound is determined from the data. 

5.1 Error bound for H w {Hl) z 1/2 Y 

To use Zolotarev approximation for the (if 2 ) -1 / 2 in the overlap Dirac operator 
(1), one needs to rescale H w to h w = H w /\ min such that the eigenvalues of h 2 w 
fall in the interval [1,6], where b = (A max /A m j n ) 2 . Explicitly 4 

1 . . d hl + C 2 l _ 1 , 2 ^ k _ /„2\-l/2 /.ox 



Amin i = i hi + C21-1 ^min ™ ; =1 ^ + 

where 

2A A 1 + C21-! 

do = I+aSt+^- (49) 

k = do ^Y'^ . (50) 



4 Here we emphasize that do (49) can be obtained explicitly, without using the condition 
: max[l - y/xr^x)}^^ = - min[l - ^/xr n (x)]\i< x < b . 
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First, we consider the multiplication of H W (H 2 ) 1 / 2 to a nonzero column 
vector Y 

(1 i n b n 
-= Y ~ h w (h 2 w + c 2n ) ]T ry-^ Y = h w (h 2 w + c 2n ) ]T 6i^i (51) 
^ H 2 J 1=1 K + C2I-1 , =1 

where the last summation can be evaluated by invoking a conjugate gradient 
process to the linear systems 

(h 2 w + ca-OZ, = y, Z = l,.-.,n. (52) 

In order to improve the accuracy of the rational approximation as well as to 
reduce the number of iterations in the conjugate gradient loop, it is advan- 
tageous to narrow the interval [1,6] by projecting out the largest and some 
low-lying eigenmodes of H^. Denoting these eigenmodes by 

H w Uj = XjUj, j = 1, • • • , k, (53) 

then we project the linear systems (52) to the complement of the vector space 
spanned by these eigenmodes 

k 

[hi + c2,_i)Z, = y = (1 - U M) Y > l = h---,n. (54) 

In the set of projected eigenvalues of H 2 , {X 2 ,j = 1, • • • , k}, we use y? max 
and y? min to denote the least upper bound and the greatest lower bound for 
the eigenvalues of H 2 , where 



k 

h w = h w — A 

Then the eigenvalues of 

h 2 = H 2 /A 2 

w wl mi 



jUjUj . 



fall into the interval (1,6), 6 = (X max /X m i n ) 2 . 

Now the matrix-vector multiplication (51) can be expressed in terms of 
the projected eigenmodes (53) plus the solution obtained from the conjugate 
gradient loop (54) in the complementary vector space, i.e., 

1 1 n k A 

H W ^=Y ~ —H w {h 2 w + c 2n ) £ b l Z l + £ ^=u 3 u)Y = S (55) 

yJH* A min l=1 j=1 y'Aj 

Then the error of S can be measured in terms of 

yty ' » ( 56 ) 
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which is zero if (55) is exact. Now assuming the errors due to the conjugate 
gradient and the projected eigenmodes are negligible compared with that due 
to the Zolotarev approximation of (hf )~ 1 ^ 2 , then we can derive the theoretical 
upper bound of er as follows. 
First, we rewrite S (55) as 

k \ N / \ \ 2 



S = J2 ~lh UjU j Y + S sign(A j ) v ^~r z (a; i )M i M}F , ^ = J (57) 



where Xj, Uj, j — k + 1, • • • , N denote the eigenvalues and eigenfunctions of H u 
Then it is straightforward to derive 



N 



S^S - yty = yt £ [(y/x]r z ( Xj )) 2 - l]u jU ) Y , (58) 



\j=k+l 



where the orthonormality ( u\uj = Sij ) has been used. Thus it follows that 



N 

,t 



\&S - Y^Y\ < k m^ <N \(^r z ( Xj )r - 1| yt ( £ uiU[ ] y 



max |(./S7r z (a; 7 -)) 2 - II y f y (59) 
where the inequality ( due to the completeness relation J2iLi u i u \ — 2 ) 

y f ( Yl u i u \ ) y < yty > (fc>o) 

has been used. Then the theoretical upper bound of (56) is 

a < max |( A /x7rz(^i)) 2 — 1| — 2 max \+/xlrz(xj) — 1\ . (60) 

fc+l<?'<JV V 1 k+l<j<N IV J y J ' ' y ' 

Thus a is less than two times of the maximum of |1 — ^Jx~jr z {xj)\ for all 
eigenvalues of h^. It follows that a must be less than two times of the maximum 
of |1 — s/xr z [x)\ for all x G [1,6], i.e., 

\&S-Y<Y\ 2(1 -A) nl , „ 
o- = J y^Y < 1 + A = z(n ' } ' (6 } 

where A is defined in (22). The inequality (61) is one of the main results of 
this paper. 

A remarkable feature of (61) is that it holds for any nonzero column vector 
y. Then one immediately sees that the deviation (6) which measures the 
chiral symmetry breaking due to Zolotarev approximation also has the same 
theoretical upper bound, 

A z <2d z (n,b) . (62) 
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Thus 2dz(n,b) not only serves as the theoretical error bound for the matrix- 
vector multiplication H W {H^)^ 2 Y ^ but also provides the upper bound for the 
chiral symmetry breaking due to Zolotarev approximation. Note that dz(n, b) 
does not depend on the lattice size explicitly. Therefore, by choosing the proper 
values of n and b ( i.e., by projecting the high and low- lying eigenmodes of 
), one can practically preserve the exact chiral symmetry of the overlap Dirac 
operator to very high precision, for any gauge configurations on a finite lattice. 

Moreover, for any gauge configuration, it is unlikely that any of the eigen- 
values of h%j would coincide with one of those 2n + 2 positions with maximum 
deviation (34), thus one usually obtains a a much smaller than the theoretical 
error bound 2dz(n, b). 

In Table 1, we list the values of a (56) for several gauge configurations 
on three different lattices respectively, along with the theoretical error bound 
2dz(n,b). It is clear that a is always much smaller than the theoretical error 
bound. Note that for each configuration, we only show the largest a among a 
set of several hundred a values, each is computed with a different Y at every 
iteration of the outer CG loop. The average value of a is usually less than 1/5 
of the largest one listed in Table 1. Details of our computation are described in 
[9]. Here we have set the stopping criterion for the inner and outer conjugate 
gradient loops at e = 10 -11 , and the error of the projected eigenmodes around 
10~ 13 . It is remarkable that a turns out to be much less than e, in contrast 
to one's naive expectation. In other words, even if the precision of the overlap 
Dirac quark propagator is only up to 10~ n , its exact chiral symmetry can 
attain A z < 1(T 12 . 

5.2 An empirical formula for the error bound 

In Table 2, we list the values of d z (n, b) = (1 — A)/(l + A) of Zolotarev optimal 
rational appproximation (29), for n ~ 10 — 20 and b ~ 10 — 10 6 . We can use 
Table 2 to determine how many Zolotarev terms ( n ) is needed in order to 
attain one's desired accuracy, after the highest and the low-lying eigenmodes 
are projected out and b = (\max/\nin) 2 has been obtained. Conversely, for a 
fixed number of Zolotarev terms, say n, one can use Table 2 to determine what 
ranges of high and low-lying eigenmodes of should be projected in order 
to attain one's desired accuracy. 

In Fig. 3, we plot d z (n, b) = max |1 — ^/xr^ ,n \x)\ = (1 — A)/(l + A) versus 
the degree n, for different values of b ranging from 10 3 — 10 6 . For any fixed 
value of b, the error bound converges exponentially with respect to n, and it 
is well fitted by 

e z (n, b) = A{b) exp{-c(6)n} , (63) 
as indicated by the solid lines shown in Fig. 3. The parameters c{b) and A{b) 
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lattice size 


(3 




^max 


b 


n 


a 


2d z {n 1 b) 


16 3 x 32 


6.0 


0.1731 


6.258 


1307.00 


12 


6.0 x ur 1 ' 2 


1.4 x 10- 1U 


16 3 x 32 


6.0 


0.1943 


6.260 


1038.01 


12 


7.0 x 10~ 12 


7.5 x 10~ n 


16 3 x 32 


6.0 


0.1767 


6.261 


1255.49 


12 


8.0 x 10~ 12 


1.2 x 10- 1U 


16 3 x 32 


6.0 


0.1955 


6.260 


1025.31 


12 


5.0 x 10~ 12 


7.2 x 10~ n 


12 3 x 24 


5.8 


0.1176 


6.210 


2788.49 


16 


7.0 x 10~ 14 


4.8 x 10~ 13 


12 3 x 24 


5.8 


0.1285 


6.211 


2336.24 


16 


6.6 x 10~ 14 


3.2 x 10~ 13 


12 3 x 24 


5.8 


0.0988 


6.206 


3945.57 


16 


2.0 x 10~ 13 


1.3 x 10~ 12 


12 3 x 24 


5.8 


0.1415 


6.213 


1927.92 


16 


5.5 x 10~ 14 


1.6 x 10~ 13 


10 3 x 24 


5.8 


0.1242 


6.214 


2503.22 


16 


1.6 x 10~ i3 


3.4 x 10~ 13 


10 3 x 24 


5.8 


0.1381 


6.209 


2021.42 


16 


6.1 x 10~ 14 


2.2 x 10~ 13 


10 3 x 24 


5.8 


0.1376 


6.204 


2032.86 


16 


5.5 x 10~ i4 


2.6 x 10~ i3 


10 3 x 24 


5.8 


0.1181 


6.204 


2759.59 


16 


4.4 x 10~ 14 


5.0 x 10~ 13 



Table 1: The error a (56) of the matrix- vector multiplication H W {H'^)'^ Y, 
for several gauge configurations on three different lattices. Evidently, each a 
is quite smaller than the corresponding theoretical error bound 2d z (n,b) = 
2(1-A)/(1 + A). 



b 


n 




10 


12 


14 


16 


18 


20 


10 


4.8 x 10~ 18 


1.9 x 10~ 21 


7.2 x 10" 25 


2.8 x 10~ 28 


1.1 x 10~ 31 


4.1 x 10" 35 


50 


1.3 x 10~ 13 


3.5 x 10" 16 


9.5 x 10" 19 


2.6 x 10~ 21 


6.9 x 10~ 24 


1.9 x 10~ 26 


100 


2.5 x 10~ 12 


1.2 x 10~ 14 


5.5 x 10~ 17 


2.6 x 10" 19 


1.2 x 10~ 21 


5.8 x 10~ 24 


500 


3.8 x 10- 1() 


4.8 x 10~ 12 


5.9 x 10~ 14 


7.3 x lO" 16 


9.0 x 10~ 18 


1.1 x 10~ 19 


1000 


2.0 x 10~ 9 


3.4 x lO" 11 


5.8 x 10~ 13 


9.8 x 10" 15 


1.7 x 10~ 16 


2.8 x 10~ 18 


2000 


8.4 x 10~ 9 


1.9 x 10~ iU 


4.2 x 10~ 12 


9.3 x 10~ 14 


2.1 x 10- Lb 


4.6 x 10" 1 ' 


3000 


1.8 x 10~ 8 


4.6 x 10" 10 


1.2 x 10" 11 


3.0 x 10~ 13 


7.7 x 10" 15 


2.0 x 10~ 16 


4000 


2.9 x 10~ 8 


8.3 x 10- 10 


2.3 x 10" 11 


6.6 x 10~ 13 


1.9 x 10~ 14 


5.3 x 10~ 16 


5000 


4.3 x 10~ 8 


1.3 x 10~ 9 


3.9 x 10" 11 


1.2 x 10~ 12 


3.6 x 10~ 14 


1.1 x 10" 15 


6000 


5.7 x 10~ 8 


1.8 x 10" 9 


5.8 x lO" 11 


1.9 x 10~ 12 


6.0 x 10~ 14 


1.9 x 10" 15 


7000 


7.2 x 10~ 8 


2.4 x 10" 9 


8.1 x 10" 11 


2.7 x 10~ 12 


9.1 x 10~ 14 


3.1 x 10" 15 


8000 


8.9 x 10~ 8 


3.1 x 10" 9 


1.1 x 10" iU 


3.7 x 10~ 12 


1.3 x 10" i3 


4.5 x 10~ Lb 


9000 


1.1 x 10~ 7 


3.8 x 10" 9 


1.4 x 10" 1U 


4.9 x 10~ 12 


1.8 x 10" 13 


6.4 x 10" 15 


1 x 10 4 


1.2 x 10"' 


4.6 x 10" 9 


1.7 x 10~ iU 


6.3 x 10~ 12 


2.3 x 10~ i3 


8.6 x lO" 115 


5 x 10 4 


9.5 x 10~ 7 


5.2 x 10~ 8 


2.9 x 10" 9 


1.6 x 10- 1() 


8.6 x 10~ 12 


4.7 x 10~ 13 


1 x 10 b 


2.0 x 10~ b 


1.3 x 10-'' 


8.0 x 10" 9 


5.0 x 10" iU 


3.2 x lO" 11 


2.0 x 10~ i2 


5 x 10 5 


8.7 x 10"° 


7.3 x 10~ 7 


6.1 x 10~ 8 


5.0 x 10" 9 


4.2 x lO" 1 " 


3.5 x 10- 11 


1 x 10 e 


1.5 x 10~ 5 


1.4 x 10- B 


1.3 x 10~ 7 


1.2 x 10~ 8 


1.1 x 10" 9 


1.0 x 10- 10 



Table 2: The deviation d z (n, b) = max ]1 - y/xr%' n) (x)\ = (1 - A)/(l + A) of 
the Zolotarev optimal rational approximation (29), versus the upper bound b 
of x G [1, b], and the degree n of the rational polynomial. 
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are determined in Figs. 4 and 5 respectively, 

c(6) = 9.17(10) ln(6)- a774(5) (64) 
A(b) = 0.465(9) ln(6) a596(9) (65) 

With the empiraical formula (63), one can estimate the theoretical upper 
bound of a or A z more conveniently, without evaluating ellitpic functions at 
all. 

5.3 Zolotarev rational polynomial of the form 7~( n-1 ' n ) 

At this point, if one recalls Chebycheff's therorem for r^ n ' m \ one may ask 
whether there exists optimal rational approximation to x' 1 ^ 2 , which has the 
form r^ n ~ 1,n \ with 1 — \fxr^ 1,n \x) having 2n + 1 alternate change of sign 
in [1,6]. Note that in the second principal n-th degree transformation (21), 
sn(x/M; A) has periods 2K/M and 2iK' /[(2n + l)M]. Thus there must exist a 
similar transformation in which sn(:r/M; A) has periods 2K/M and 2iK'/2nM. 
Explicitly, it reads 



\m ) v m n#=i 1 + xlu: k /C 2 /_i 



where 



mnr=i[i + ^(«;«)/C2j-i] 



2n Q2(m£L. K > 

m = ' ( ' 

l_ sn 2(^ ;K A 



V 2n ' 

Ci = ; ^f^'L (69) 

2n ' 



Then the Zolotarev optimal rational polynomial of the form r*™ 1,n ) is 

Jn-Ln),, _ 2A 1 Y[^{l+X/C 2l ) 

z W_ i+Amnr =1 (i+^ 2M ) UUj 

The difference 1 — y/xr^~ 1,n \x) has 2n + 1 alternate change of sign in [1,6], 
( b = k~ 2 ) , and it attains its maxima and minima alternatively as 

1- A 1- A 1 - A 1 - A 



1 + A' 1 + A' ' 1 + A'l + A 
at 2n + 1 successive points x E [1,6] : 

1 1 



1, 



1 _ kW [K_. K ,y ' ! _ «,2 sn 2 «') ' 

14 



K 2 



Obviously, for any given n and b, the deviation of 1 ' n " > is larger than 



that of r^' n \ i.e., 



d{r [ r hn) ) = l YVK > W\ = d{r(r)) ' (71) 
and for most cases, 

d{r { r 1,n) ) ^ 2.5 x d(rp n) ) (72) 

In Table 3, we list the deviation d(4™ _1 ' n) ) = (1 - A)/(l + A) of the 
Zolotarev optimal rational appproximation (70), versus the degree n, and the 
upper bound b of x G [1,6]. Comparing the corresponding entries ( with the 
same b and n ) in Table 3 and Table 2, one immediately sees that d{r%~ 1,n) ) ^ 
2.5 x d(rp n) ). 

If one uses (70) to approximate the inverse square root of in the overlap 
Dirac operator, then one has 

1 _ A. n?"! 1 ^ + C7 a ) D ^ B t (?3) 



where 

D ° ~ i + A nr^(i + ^) (74) 

= A ?T P^^% j , • (75) 

Comparing (48) with (73), one immediately sees that it is more advanta- 
geous to use the former approximation than the latter, especially for computing 
quark propagators, since only one more matrix multiplication with (h^ + c 2n ) 
after the completion of the inner CG loop would yield about 2.5 times higher 
accuracy than using (73). Although we used (73) for our computations in Ref. 
[9], we have switched to (48) for better accuracy, in our ongoing lattice QCD 
computations. 



6 Concluding remarks 

In this paper, we have discussed the basic principles underlying the rational 
approximation, and shown explicitly that the Zolotarev approximation is in- 
deed the optimal rational approximation for the inverse square root function. 
For the overlap Dirac operator, we have derived the theoretical error bound 
for the matrix- vector multiplication H^H^^^Y ', which is equal to two times 
of the maximum deviation dz(n, b) of the Zolotarev rational polynomial. This 
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b 


n 




10 


12 


14 


16 


18 


20 


i nnn 

1UUU 


0.0 A 1U 


Q A v 1 n -11 


1.0 A 1U 


9 7 v 1 n~ 14 

Z. / A 1U 


^.0 A 1U 


/ .O A 1U 


9nnn 

zuuu 


2.2 x 10~ 8 


4.8 x 10" iU 


1.1 x lO^ 11 


2.4 x 10~ i3 


5.3 x 10~ ib 


1.2 x 10~ lb 


3000 


4.5 x 10~ 8 


1.1 x 10~ 9 


2.9 x lO" 11 


7.5 x 10~ 13 


1.9 x 10~ 14 


5.0 x 10~ 16 


4000 


7.2 x 10~ 8 


2.0 x 10~ y 


5.7 x lO" 11 


1.6 x 10" 12 


4.6 x 10~ 14 


1.3 x 10~ lb 


5000 


1.0 x 10~ 7 


3.1 x 10~ y 


9.4 x 10~ n 


2.8 x 10~ 12 


8.6 x 10~ 14 


2.6 x 10~ 15 


6000 


1.3 x 10~ 7 


4.3 x 10~ y 


1.4 x 10- 10 


4.4 x 10~ 12 


1.4 x 10~ 13 


4.5 x 10~ 15 


7000 


1.7 x 10~ 7 


5.7 x 10~ y 


1.9 x 10" 10 


6.4 x 10~ 12 


2.1 x 10~ 13 


7.2 x 10~ 15 


8000 


2.1 x 10~ y 


7.1 x 10~ y 


2.5 x 10~ iU 


8.7 x 10~ i2 


3.0 x 10~ 13 


1.1 x 10~ 14 


9000 


2.4 x 10~ 7 


8.7 x 10~ y 


3.1 x 10" 10 


1.1 x lO^ 11 


4.1 x 10~ 13 


1.5 x 10~ 14 


10000 


2.8 x 10~ 7 


1.0 x 10~ 8 


3.9 x lO" 10 


1.4 x 10" 11 


5.3 x 10~ 13 


2.0 x 10~ 14 



Table 3: The deviation <i(r^'~ 1 ' n ' ) ) = max |1 — y^xr^ l,n \x)\ = (1 — A)/(l + A) 
of the Zolotarev optimal rational approximation (70) versus the upper bound 
b of x e [1, b], and the degree n of the rational polynomial. 

is also the upper bound for the chiral symmetry breaking due to Zolotarev 
approximation. Some numerical values of dz{n, b) are listed in Table 2 as well 
as plotted in Figure 3. An empirical formula for dz(n, b) is determined, which 
provides a reliable estimate of the theoretical error bound, especially for the 
range of parameters : b ~ 10 3 — 10 4 and n ~ 10 — 20. We also compare 
two possible forms of Zolotarev optimal rational approximation : (x) and 
r^z~ 1,n \x), and point out that the former seems to be the better choice for 
computing quark propagators, since with the same computational cost, one 
has d(r^ ,n ^) — 0.4 x rf(r| l ~ 1 ' n ' ) ). 

With Zolotarev optimal rational approximation for (iJ^)~ 1//2 in the overlap 
Dirac operator, one has no difficulties to preserve exact chiral symmetry to very 
high precision ( e.g., A z < 10~ 12 ), for any gauge configurations on a finite 
lattice. This feature is vital for lattice QCD to extract physical observables 
from the first principles. In practice, one might have difficulties to push the 
error o (56) down below 10~ 13 , which is essentially due to the inaccuracies 
of the projected ( high and low-lying ) eigenmodes of H 2 , rather than the 
Zolotarev approximation of (H 2 )^ 1 ^ 2 . In the future, we will try to improve 
the accuracy of the projected eigenmodes. In the meantime, the precision of 
exact chiral symmetry up to 10~ 13 should be sufficient for many calculations 
in lattice QCD with overlap Dirac quarks. 

This work was supported in part by the National Science Council, ROC, 
under the grant number NSC90-2112-M002-021, and also in part by NCTS. 
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Figure 1: The difference 1 — y/xrz(x) of Zolotarev optimal rational approx- 
imation (29) with n = 6 and b = kT 2 = 1000. Note that there are exactly 
2n + 2 = 14 alternate change of sign in [1, 1000], and the maxima and minima 
have exactly the same magnitude. 
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Figure 2: The positions Xi (34) of the maxima and minima of 1 — yfxrz(x) for 
n = 20 and b = k~ 2 = 6000. 
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Figure 3: The error bound of Zolotarev optimal rational approximation for the 
inverse square root function x" 1 / 2 , x G [1, b], versus the degree n of the rational 
polynomial, and for several different values of b ranging from 10 3 to 10 5 . For 
any fixed value of b, the error bound converges exponentially with respect to 
n, and is well fitted by A{b) exp{-c(6)n} with c(b) = 9.17(10)ln(6)-°- 774 ( 5 ) and 
A(b) = 0.465(9)ln(6) a596 ( 9 ). 
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Figure 4: Determination of the parameters of c(b) 
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Figure 5: Determination of the parameters of A(b) 
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